Light propagation in finite-sized photonic crystals: 
Multiple scattering using an electric field integral equation 
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We present an accurate, stable and efficient solution to the Lippmann-Schwinger equation for 
electromagnetic scattering in two dimensions. The method is well suited for multiple scattering 
problems and may be applied to problems with scatterers of arbitrary shape or non- homogenous 
background materials. We illustrate the method by calculating light emission from a line source in 
a finite sized photonic crystal waveguide. 



I. INTRODUCTION 

Design and development of future optical components 
rely heavily on adequate descriptions of light propagation 
in micro and nano structured environments. Although 
the governing differential equations have been known for 
over a century, analytical solutions are available only for 
a limited number of highly symmetrical problems. There- 
fore, in most cases of practical interest, numerical meth- 
ods are employed. Ideally, the numerical method should 
be accurate enough to capture all relevant physics and al- 
low for understanding the different scattering channels; 
yet it should be fast enough to ensure acceptable run- 
times for practically relevant structures. 

The invention of photonic crystals [l], 0, 0] has of- 
fered hitherto unprecedented control of light propaga- 
tion. Consequently, these novel materials are of great 
importance in the design of optical components for fu- 
ture information technology as well as fundamental solid 
state quantum optics experiments. In the context of light 
propagation in micro and nano-structured media, such 
as photonic crystals, the full wave nature of the electro- 
magnetic field needs to be taken into account. A large 
number of different methods are being explored in the 
investigation of light propagation in these materials, in- 
cluding plane wave expansion for band structure calcula- 
tions [j] , the Generalized Multipole Technique [f| and the 
Finite Difference Time Domain (FDTD) method [f|. In 
particular, we note that for calculations involving large 
numbers of scatterers in an otherwise homogeneous back- 
ground, Rayleigh-multipole methods have been used for 
calculations on microstructured fibres @, H, 
as photonic crystals composed of cylinders 
spheres (T5 , The use of multipole expansions for 

the fields ensures a significant reduction in the number 
of basis functions, thus enabling calculations on complex 
structures of practical interest. 

In addition to the above methods, various volume and 
surface integral methods exist for the electric and/or 
magnetic fields. These methods typically rely on dis- 
cretization to express the integrals as linear systems of 
equations - a procedure known as the method of mo- 
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ments. In the discretization process the solution is 
expanded in terms of linearly independent basis func- 
tions which vary in complexity, depending on the specific 
method employed [l5j]. A particularly simple and conve- 
nient choice of basis function is the pulse basis function 
resulting in what is usually termed the coupled (or dis- 
crete) dipole approximation (CDA) [H, [Oj- The CDA 
allows for relatively easy implementation as well as the 
physically attractive property that the resulting field can 
be directly interpreted as the sum of the field from all 
the individual cells in the scattering structure oscillating 
as dipoles in response to the incoming field. Depending 
on the desired accuracy and the nature of the scatter- 
ing problem at hand, however, the required number of 
basis functions may become too large for practical calcu- 
lations on e.g. photonic crystals. In addition, this type of 
discretization may lead to stability problems which, for 
example, limit the efficiency of volume integral equations 
for the electric field in the case of high refractive index 
contrasts [i~5j . 

For use in the modelling of nanophotonic structures 
and in particular in the context of light-matter interac- 
tion in photonic crystals, we focus in this paper on an 
integral type scattering formulation of light propagation 
based on the electric field Green's tensor. The Green's 
tensor, G(r, r'), is the electromagnetic propagator which, 
for a given system, holds all information necessary to 
solve the inhomogeneous vector Helmholtz equation and 
may be interpreted as the electric field at point r due 
to an oscillating point dipole at the point r'. Of special 
importance in nanophotonics modelling is the imaginary 
part of the Green's tensor at r = r' which is propor- 
tional to the local optical density-of-states (LDOS) [18l |. 
The LDOS and/or the Green's tensor for photonic crys- 
tals have previously been investigated theoretically us- 
ing plane wave expansions [H Hi and FDTD [H [H 
as well as through an expansion in eigenstates [23l . |24| . 
The application of the CDA to the calculation of the 
Green's tensor in a photonic crystal slab was reported in 
Ref. [25|. Using the Green's tensor for the background 
medium, one may recast the wave equation as a scatter- 
ing problem, in which case the solution is given in terms 



2 



of the so-called Lippmann-Schwinger integral equation 
(26| . Since the Green's tensor is known for a number of 
simple material configurations, these may be used as the 
background material in order to limit the extend of the 
numerical calculations. In particular, expressions for the 
Green's tensor for stratified media are available (27l . [28| . 

In this article we describe a multiple scattering solu- 
tion to the Lippmann-Schwinger equation. The method 
is developed in two dimensions, but we note that a sim- 
ilar approach is possible for three dimensional problems 
as well. The method will find applications in modelling 
of nanophotonic structures and devices, such as waveg- 
uides, junctions and filters as well as switches and single 
photon sources based on photonic crystals. In light of 
the above discussion the method may be viewed as a 
hybrid between integral type method of moments calcu- 
lations and multiple scattering multipole methods. In 
our approach, the Lippmann-Schwinger equation is first 
expanded in cylinder functions (so-called normal modes) 
and solved within the scattering elements. Solutions at 
points outside the scattering elements are subsequently 
calculated directly using the Lippmann-Schwinger equa- 
tion. Because of the integral formulation the method 
may benefit from known results for the Green's tensor 
in the background material while the normal mode ex- 
pansion reduces the number of basis functions needed, 
thus enabling calculations on material systems of prac- 
tical relevance. In addition, we make use of a number 
of theorems which are regularly employed in multipole 
methods in order to simplify the evaluation of the scat- 
tering matrix elements. 

The article is organized as follows: Section |TT] intro- 
duces the method and provides a number of results that 
may be of value in the practical implementation. The 
calculations are illustrated in section lTTFl through the so- 
lution of two small scattering problems and in section |TTT] 
we use these as examples to discuss a practical method of 
evaluating the accuracy of a given calculation. In section 
IIVI we provide an example application of the method in 
the form of a two-dimensional photonic crystal waveg- 
uide at the edge of a dielectric block. Section IVl gives the 
conclusions. 



II. FORMULATION OF SCATTERING 
PROBLEM AND SOLUTION METHOD 

We consider scattering of monochromatic light, 
E(r, t) — E(r) exp(— iojt) in two dimensions, correspond- 
ing to scattering problems in which the geometry is in- 
variant along the ^-direction and the light travels in the 
xy-plane. We limit the discussion to non-magnetic ma- 
terials and consider general material systems, consisting 
of a finite number of piecewise constant dielectric pertur- 
bations to a piecewise homogeneous background. This is 
the case, for example, for most photonic crystals. The 
electric field E(r) solves the vector wave equation 

V x V x E(r) - fc^e(r)E(r) = 0, (1) 



where e(r) is the position dependent permittivity and 
fco = |ko| = uj/c is the wave number in vacuum. For two- 
dimensional problems, the vector equation decouples into 
two independent equations for the Transverse Electric 
(TE) and the Transverse Magnetic (TM) polarizations. 
In the case of TE polarization, the electric field is ori- 
ented in the xy-plane, whereas, for TM polarization the 
electric field is oriented along the z-axis and the scatter- 
ing calculation is essentially a scalar problem. In order 
to reformulate Eq. ^ as a scattering problem, we con- 
sider the change in permittivity, Ae(r) = e(r) — £s(i"), 
caused by the introduction of scattering sites into the 
background medium described by cb (r) . The solution to 
Eq. (P) with e(r) = e#(r) is denoted by E s (r) and rep- 
resents the incoming field. The full solution to Eq. JTJ) is 
the sum of the incoming field and the scattered field. It 
is given by the Lippmann-Schwinger equation 

E(r) = E B (r)+ / G s (r, r') fcg Ae(r') E(r')dr', (2) 
Jd 

in which G s (r, r') is the Green's tensor for the back- 
ground medium. For a homogeneous background, 
es(r) = es = n B , the two dimensional Green's tensor 
is given as [291 ] : 

G^(r, r') = (l + {H (k B \r r'|) (3) 

in which ks = «sfco is the wave number in the back- 
ground medium and Ho(r) is the Hankcl function of the 
first kind of order 0. 



A. Expansion in normal modes 

In two and three dimensions the real part of the 
Green's tensor diverges in the limit r = r', which means, 
that for integrals in which r is inside the scattering vol- 
ume (such as in this work) an alternative formulation of 
the Lippmann-Schwinger equation must be employed in 
which the singularity is isolated in an infinitesimal prin- 
ciple volume and treated analytically [3(|. Therefore, we 
follow Ref. [29[ and rewrite the equation as 

E(r) =E B (r)+ lim / G B (r, r') fcg Ae(r')E(r')dr' 
-L^E(r), (4) 

(B 

in which L xx = L yy = 1/2 and L a p = otherwise, corre- 
sponding to a circular exclusion volume 6V centered on 
r' = r. 

We will solve for the total field, E(r), inside the scatter- 
ing material only, as the solution everywhere else can be 
subsequently obtained by use of the Lippmann-Schwinger 
equation. The incoming field, E B (r), is a solution to 
the wave equation with no scatterers and thus, in gen- 
eral, can be expanded on the solutions in the bulk back- 
ground material. Similarly, the total field at positions 
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inside each scatterer may be expanded on the solutions 
in a homogeneous material with the same permittivity. 
Based on these considerations we construct a basis con- 
sisting of normal modes, each with support on only one 
of the scattering sites: 

= K n J q (k d r d y^ S d (r) 
^ = K*J q (k B r d )e l «^S d (r), 

where we have defined a combined index n = (q, d, a), in 
which q denotes the order of the normal mode and d de- 
notes the particular subdomain D d of the scattering ma- 
terial we consider, e.g. which cylinder. For convenience 
we have included the field component a G {x, y, z} in the 
index n as well. J q {r) denotes the Bessel function of the 
first kind of order q and k d — n d ko, where n d denotes the 
refractive index in the subdomain D d . (r d ,ip d ) denote 
polar coordinates with respect to the local coordinate 
system in the subdomain D d and S d = 1 for r G D d and 
S d = otherwise. K n and are normalization con- 
stants. To ease the notation we will write n in place of 
any of the indices q, d, a and we will supress the index on 
the coordinates (r, tp) as this will not lead to ambiguities. 
We define an inner product as 



(ipm\il>n) = / ft(r)^n(r)dr, 



(5) 



which is in general non-zero for m ^ n, and we normalize 
the basis functions so that {ip n \ip n ) = (ipnl^Pn) = 1- By 
expanding the electric fields as 

E(r) = y^e n V> n (r)e n 

n 

E» = £e^(r)e„, 

n 

where e„ is a unit vector in the direction a, and pro- 
jecting onto the basis formed by ip m e m and ip^e m , the 
Lippmann-Schwinger equation may be rewritten in ma- 
trix form as 

A< 



^ n n 




P -D 



•' Rd 



(6) 



FIG. 1: Sketch of the local coordinates used for the calculation 
of the self-term in scattering domain D. 



is very sparse, since basis functions belonging to different 
scattering domains are orthogonal by construction). The 
underlying strategy and the practical implementation of 
the two methods, however, are very different. In the 
CDA, the projection onto the basis functions is straight 
forward, resulting typically in a very large matrix equa- 
tion system that is subsequently solved by iterative meth- 
ods. In the present method it is the projections onto the 
basis functions that are potentially time consuming, but 
eventually leads to a relatively small system of equations. 
A number of mathematical results exist, though, that can 
dramatically speed up the calculations of the matrix el- 
ements and the implementation in general. In sections 
III BIIII El we discuss these results. 

The matrix elements need to be evaluated for basis 
functions within the same domain (self-terms) as well as 
different domains (scattering terms). The case of a ho- 
mogeneous background is of special interest as one will 
often be able to separate the background Green's tensor 
into terms corresponding to a homogeneous background 
and a number of additional scattering terms. Therefore, 
we focus in this section on the evaluation of the matrix 
elements for a homogeneous background Green's tensor, 
Eq. ([3]), and return to the additional terms due to scat- 
tering in an inhomogeneous background in section [TV] 



Self-terms 



in which 

Gtfn= [ C>(r) km / G^(r, r>„(r')dr'dr, 

JD dV->0 J D _ sv 

(7) 

and we have written explicitly the directions a, corre- 
sponding to the indices m, n for clarity. 

Eqs. © - © hold the complete formulation of the 
method. The form of the central matrix equation © is 
only slightly different from that of the CDA in that the 
left hand side is generally a matrix (although this matrix 



The evaluation of the self-terms is complicated by the 
divergence in the Green's tensor for r' = r and the fact 
that the integrand couples r' and r, effectively resulting 
in a four-dimensional integral. In the following we adopt 
a particular method for the evaluation, in which the four- 
dimensional integral, Eq. ©, is rewritten in terms of a 
number of one and two dimensional integrals. Although 
this method is suitable for evaluation of general matrix 
elements, we note that other methods may be more effi- 
cient in case of a particular geometry and/or polarization. 

Figure Q] shows a sketch of the local coordinates used 
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for the evaluation of the integral. In order to efficiently 
treat the divergence, for each r we first integrate r' over 
the entire plane P less the principal volume centered on 
r. Subsequently, we subtract the integral for r' ^ D, for 
which the limit SV — > is trivial, since r' ^ r. The 
matrix element is thus rewritten as G^ n — A^ n — B^f n , 
in which 

A a ±= I r m (r) lim / Gf^r^OVnCrOdr'dr 

JD dV^O j p_ SV 

and 

B«i = / / G^(r,r')Vn(r')dr'dr. 

Using the Graf addition theorem (cf . appendix [X| , 
we can simplify the evaluation of A a ^ by expanding the 
function ip n (r') around the point r, so that 

A af} = J2 [ J m (k R r)J n+ ^(k R r)e^- m ^ 

x / Gf /3 (r,R)J, I (fc i? i?)(-l)^e-^dRdr 
Jp-sv 

= V / J m (fc fl r)J„ + , I (fc i? r)e l (" + ^ m ^dr x if, 

where k R — ridko and we have exploited the fact that 
the integration over R is over the entire plane (less the 
principal volume), and thus does not depend on r. The 
evaluation of this integral over R, denoted by if above, 
can be reduced to a number of one dimensional integrals 
as shown in appendix[Bj We note that the simple angular 
dependence of the integrands in many cases allows for a 
reduction of the remaining integral over r to a sum of 
one dimensional integrals. 

Evaluation of B%f n may also be substantially simplified 
using the Graf addition theorem to expand the Hankcl 
function in terms of Bessel and Hankel functions defined 
in the local coordinate system. The expansion differs 
depending on the sign of r — r'; for r' > r we write the 
integrand as 

6^ n (r ) r')=C(r)G^(r,r')^„(r / ) 

= Jm(k R r)e' lm ^C^ l l -J^k B r)e-^\ 

x H^k B r , )J n (kRr'e l{n+ ^'), (8) 
whereas, for r' < r we write 

&^(r,r') =J2Jm(k R r)e- im ^£^ { l -H,{k B r)e^\ 

x J^kBr'^kRr'e 1 ^^'), (9) 

in which C a @ is the a, (3 component of the linear op- 
erator in Eq. Derivatives for the general cylinder 



functions d a dp{Z fl (k B r)e ± ' llJ,v '} are provided in appendix 
[Q For circular scatterers we always have r 1 > r and 
the expression for 23£m factors into a number of one di- 
mensional integrals. Similarly, the evaluation of B^ n for 
non-circular scatterers may be conveniently split depend- 
ing on whether r' is outside or inside the circumscribing 
circle (denoted by Ru in Fig. [1}. In the former case, 
the expression factors into separate integrals for r and 
r', whereas in the latter case, the two integrations are 
coupled. Again, the simple angular dependence of the 
integrands in many cases allows for a reduction of these 
integrals to a sum of two dimensional integrals. 



C. Scattering terms 




FIG. 2: Sketch of local coordinates for r and r' in two in- 
dependent scatterers. Using the Graf addition theorem the 
calculation of the scattering matrix elements may be simpli- 
fied considerably and expressed in terms of the integral over 
the individual scatterers and the distance A between the cen- 
ters of the two local coordinate systems. 

For the calculation of scattering terms the integration 
domains for r and r' are completely separated in space 
and so the Green's tensor is well behaved at all points of 
interest. In this case we employ the Graf addition theo- 
rem twice to express the Hankel function in terms of the 
distance between the centers of the two local coordinate 
systems, as illustrated in Fig. [2] 




x C a>3 {j x {k B r)e- lXv }Av 

x / J^k B r').J n (k n r'y( n -^*'Ar'. (10) 

Eq. (|10p shows that the scattering matrix calculation 
factors into terms that depend only on the geometries of 
the individual scatterers and the distance between them. 
Since the Hankcl function as well as the Bessel functions 
are well behaved at all points of interest, the integrals 
may be directly evaluated for arbitrary geometries, e.g. 
using numerical quadrature. Note that the procedure 
outlined above is compromised when A < R m +R n , where 
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R m and R n are the radii of the circumscribing circles of 
domains D m and D ni respectively. This could happen in 
the case of close non-circular scatterers. In this case the 
Graf addition theorem is not valid and one will need to 
employ a strategy similar to Eqs. ([8]) and ([9]). 



D. Background electric field 

The incident background electric field, E s (r) is a so- 
lution to the wave equation without the scatterers. In 
the case of a bulk background, the solutions are plane 
waves, and the expansion in terms of cylinder functions 
is readily obtained using the Jacobi-Anger identity, as 
discussed in appendix[Al Instead of using plane waves as 
the background electric field we may use the columns of 
the 2D Green's tensor. These are related to the electric 
field at r due to a line source at r' (2i|. By comparing 
with Dyson's equation 

G(r,r') =G B (r,r') 

+ / G B (r,r")fc 2 Ae(r")G(r",r')dr", 
Jv 

we see that the solution to the Lippmann-Schwingcr 
equation in this case exactly produces the corresponding 
columns of the full 2D Green's tensor for the scattering 
problem. 



E. Exterior solution 

The matrix equation (J5]) is solved using standard linear 
algebra routines to yield the solution at any point inside 
the scattering domains. The solution at any point out- 
side the scatterers can be subsequently obtained directly 
from the Lippmann-Schwinger equation which is now an 
explicit equation: 

E(r)=E»+ f G B (r,r')fc 2 Ae(r')^e„Vr l (r')e„dr'. 

(11) 

The sum in Eq. runs over all basis functions in 

all scattering domains. Again, the calculation may be 
simplified considerably by the use of the Graf addition 
theorem to rewrite the equation in terms of the distances 
from the centers of the local coordinate systems. Con- 
sidering for simplicity the case of just a single scattering 
domain D we rewrite the equation as 

E(r) —E B (r) + l -Ae d k 2 n £ C {H^k B L)e^ 8 } (-l)"e„ 
x / J AI (fc s r')e„if„J n (fc fl r')e I ("-^ v 'dr', (12) 

JD 

where now (L, 0) are polar coordinates of O' with respect 
to r, cf. appendix IA1 



F. Example calculations 

To illustrate the method we consider now an example 
scattering problem in which a TE plane wave is incident 
from the top left on a small crystallite of air cylinders 
in a high-index dielectric. Fig. [3] shows the absolute 
square of the total field as a function of the position in 
the rry-plane. Also we show the magnitude of the E x 
and the E y components of the field along the line y = 
through the centres of three of the cylinders. Clearly, the 
x component shows a number of discontinuous jumps, 
whereas the y component is continuous in accordance 
with the boundary conditions. We note that the multiple 
scattering from the air cylinders act to partly block the 
light, resulting in the formation of a standing wave in 
the upper left part of the figure. Typically we use the 
same number of basis functions in each scattering domain 
and for each polarization, so that |g| < M max . This 
calculation was performed using M max = 10, resulting in 
a matrix equation system of 294 unknowns. 

As noted in section !!! Dl we may use the present method 
to calculate the Green's tensor for the scattering struc- 




FIG. 3: Example calculation. A TE plane wave of unit am- 
plitude, E s (r) = eexp(in.gko • r), is incident from the top 
left on a crystallite consisting of seven air holes (n<j = 1) in 
a high-index dielectric background (tib = 3.5). Parameters 
are ko = (\/3/2, —1/2) and R = 0.3a where R is the radius of 
the cylindrical holes and a = 0.3Ao is the distance in between. 
Top: Absolute square, |E(r)| 2 , of the resulting field as a func- 
tion of position in the xy plane. Bottom: Absolute value of 
the components E x (x) (red solid line) and E y (x) (blue dashed 
line) along the line y = 0. 
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FIG. 4: Real (top) and imaginary (bottom) part of the to- 
tal TM Green's tensor G 22 (r, r') as a function of r with 
fcor' = (—1,1/4) (as indicated by the red dot) in a struc- 
ture consisting of four dielectric rods (rid = 3.5) of square 
cross section in air. Parameters are a = 2L where a is the 
distance between the rods and L — Ao/4 is the side length. 

ture. In Fig. [3] we consider a geometry consisting of four 
square dielectric rods in air and we show the real and 
imaginary part of the TM Green's tensor G zz (r,r') as a 
function of r for constant r' indicated in the figure. The 
real part diverges in the limit r — ► r', whereas the imag- 
inary part is continuous at all points. In the limit r = r' 
it is proportional to the LDOS, as noted in the introduc- 
tion. The calculaton was performed using M max = 10, 
resulting in only 44 unknowns. 



where the integrals are taken over the area of the scatter- 
ing sites only. Fig. [5] shows the global error as a function 
of the number of basis functions used in the expansions 
and dependent on the error in the matrix elements for the 
solutions depicted in Figs. [3] and [4] The error analysis 
was performed by first calculating the matrix elements to 
high precision, using an absolute error tolerance on the 
numerical integrals of 10~ 6 . Subsequently, for each value 
of M max the corresponding linear equation system, Eq. 
(|rj|) was constructed and a random complex number of 
fixed modulus, SG mn was added to each element in the 
matrix of modulus larger than SG mn before solving the 
equation system. 

The analysis shows an exponential like decrease in the 
global error as a function of the number of basis func- 
tions, underscoring the massive reduction in basis func- 
tions due to the expansion in normal modes when com- 
pared to conventional discretization methods. This is the 
case for the cylindrical holes in Fig. [3] as well as for the 
square rods in Fig. 2J The convergence is faster in the 
case of the cylindrical holes, which is partly because the 
basis functions have the same symmetry as the scatterers 
and partly because the plane wave field is easier to ap- 
proximate than the (divergent) Hankel function. Clearly, 
the artificial error on the matrix elements acts to limit 
the minimum achievable global error, and the analysis 
thus confirms that the global error is controlled by the 
number of basis functions as well as the accuracy of the 
numerical quadrature. We note, that the measure, Eq. 
(TT3")) . may be viewed as a test of self-consistency of the 
method which is of principal importance for any solution 
to Eq. @. From Fig. [5] we can see that the measure is 
also of practical importance, since, for a given tolerance 
on the numerical integrals, it can be used to estimate the 
number of basis functions needed to reach the minimum 



III. ACCURACY OF THE METHOD 

The numerical error stems primarily from evaluation of 
the matrix elements and the truncation of the basis set. 
After solving the linear equation system, Eq. ©, the 
accuracy of a given solution may be estimated by substi- 
tution back into the Lippmann-Schwinger equation. To 
this end we define the local error as 



£ L (r) = |E s (r)-E„ um (r) + / G s (r, r') fc 2 Ae(r') E num (r' 



(13) 

and we note that since E B (r) and G B (r,r') are known 
analytically, we can use this as a measure of the accuracy 
of a given solution even if we do not know the analytical 
solution. Based on the local error, we define the global 
relative error as 




JgL(r)dr 
J|E s (r)|dr : 



FIG. 5: Global error as a function of the number of basis func- 
tions used in the expansion of the electric fields (controlled by 
M m ax)- Circular markers correspond to the problem in Fig. 
13 with different curves corresponding to different fixed errors 
on the relevant matrix elements as indicated. Square markers 
corresond to the problem in Fig. [4] calculated for the Green's 
tensor (G 22 ) and plane waves (PW) as the background field. 
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global relative error. 



IV. EXAMPLE APPLICATION: LIGHT 
EMISSION IN A FINITE SIZED PHOTONIC 
CRYSTAL WAVEGUIDE 

As an example of the utility of the method we present 
in this section results for the investigation of light propa- 
gation near the edge of a finite sized two dimensional pho- 
tonic crystal. We consider a photonic crystal made from 
80 circular rods of refractive index n c i — 3.4 in a lower- 
index background (ub = 1.5). The cylinders are placed 
in a square lattice, and a short waveguide is created by 
the omission of 4 rods along the (11) direction of the crys- 
tal. The waveguide along with the crystal is terminated 
by an interface to air. We focus on TM polarized light 
and calculate the Green's tensor of the system G zz (r, r'). 
Although the integral expressions become larger, a simi- 
lar procedure as the one outlined below may be used for 
the calculation of TE polarized light as well as for mul- 
tiple interfaces. We start by extending the formalism of 
the previous sections to the case of a non-homogeneous 
background Green's tensor in section IIV Al and go on to 
show example calculations of light emission from a finite 
sized photonic crystal in section HV Bl 



A. Additional scattering near interface 

For the scattering calculations near a dielectric inter- 
face we use the Green's tensor for the dielectric half-space 
as the background Green's tensor in Eq. ([2]). The 2D 
Green's tensor for general strattified media is given in 
Ref. [28| . It is expressed in terms of an integral in k-space 
and below we discuss the calculation of the elements G^ n 
in the special case of a single dielectric interface. We con- 
sider TM polarized light incident on an interface at y = 
between two media with refractive indices ua and ns, re- 
spectively. We will deal only with scatterers in the lower 
layer (layer B) and so, following Ref. [2S| . the 2D Green's 
tensor is given as 



Gf z (r,r') = -§<5(R) 



/OO -i 
-oo k B,y 



oo ^a,y 
oo -pS 



BA c i k x (x-x') c -i k B .y(v+v') jfc 
47T k B ,y 



(14) 



where ki y = y/ kf — k x with ki = n/fco, I £ {^4, B} and 



T b - 

S BA ~ 



k B,y k A,y \^B V k A ^£ 



kB,y + k A , y ^Jk 2 B -kl + y/k A - kl 
is the Fresnel reflection coefficient. 



In Eq. (fT4f the first two terms correspond to the 
Green's tensor of the homogeneous material whereas the 
last term gives the reflection off the interface. This means 
that the evaluation of the matrix element Gi^ n naturally 
splits into a direct homogeneous material part and an 
indirect interface scattering part. The former is exactly 
what was handled in sections III Bl and III CI so we concen- 
trate in this section only on the scattering contribution 
qS . 



<~ T r, 



J~B A {k x ) 



4tt 7_ 00 k B , y (k x ) j Dm j D 



x e 



^* n {r)e lk * {x - x,) 

) JD m JD n 

<fcs.y(*.)(»+w')^ n ( r ')dr'drdfe. (15) 



In order to carry out the integration we first write 
(x,y) = (X, Y) + (r cos r sin <p) and (x',y') = (X',Y')+ 
(r' cosip' ,r' siny/), where (X, Y) and (X',Y') denote ab- 
solute coordinates of the centers of the local coordinate 
systems. We then recast the expression in terms of local 
coordinates as 

r s _ J_ f°° F BA (k x ) i(k x (x-x')-k B .vQr+Y')) 
x / J m (k m r)e-" nv e lkBr cos{v ~ e) dr 

JD m 

x / J n (k n r')e lnv ' 'e lkBr ' ' cos ^'~ e ^dr'dk x , (16) 

where we have rewritten the inner products of the wave 
vectors and the position vectors in the two domains in 
terms of the angles between them. This angle becomes 
imaginary whenever k 2 > k B . As in the case of the 
homogeneous background we are able to simplify the ex- 
pression further by factoring out the integrals over the 
domains D m and D n . To this end we use the Jacobi 
Anger identity, cf. appendix [XI to rewrite the matrix 
elements as 



mn 47,- 



;A+7 



A,7 
oo rS 



F BA {k x ) i(kJX-X')-k B . v (Y 



l kB,y{k x 

x / J m (k m r)Jx(k B r)e^ x - m ^dr 



(k ]1 {x-x')-k Bty (Y+Y')) e -i(\e+ 1 e') Ak 



x / J n (k n r')J x (k B r')e l{ ~t +n ^' dr'. 



(17) 



Due to the circular symmetry, the angular integrations 
over the domains D m and D n lead to non-zero values only 
for A = m and 7 = n. In these cases the radial integrals 
have well known analytical values, leaving only a final 
integration over k x . 
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FIG. 6: Top: Absolute value \G zz (r,r')\ of the TM Green's 
tensor for a finite sized photonic crystal waveguide consisting 
of 80 rods of refractive index = 3.4 in a background with 
an interface between a low-index dielectric (jib = 1.5) and 
air (tia = 1). The results are calculated as function of r 
with kor' = (0, —7.58) (indicated by the red dot and vertical 
dashed line). Bottom: Real (red) and imaginary (blue) parts 
of G zz (y,r') along the line x = 0. 



B. Light emission in finite sized photonic crystal 
waveguide 

In Fig. [6] we show a contour plot of the absolute 
value of the Green's tensor \G zz (r, r')| along with real 
and imaginary parts at positions along the x-axis. Re- 
sults are shown for kor' = (0, —7.58), in the center of 
the waveguide at the location of one of the missing rods. 
In an infinite waveguide, this would be the location of 
the field antinode of the waveguide mode. The periodic 
Bloch-mode character of the waveguide mode is evident 
also in the case of this finite waveguide and the structure 
acts as a resonator, greatly increasing the absolute value 
of the Green's tensor for positions r inside the waveg- 
uide as compared to the bulk medium. For r — * r' the 
real part of G zz {r, r') diverges. This is the case also in 
Fig. [6l but the divergence is too weak to show up at 
the chosen discretization. Although the finite waveguide 
acts as a resonator, light can propagate out of the end 
facet. Fig. [7] shows \G zz (r,r')\ at positions outside the 
structure. As noted in section Hi Dl the Green's tensor is 
related to the electric field at point r due to a line source 
at point r'. Therefore we may interpret the figure as the 
emission pattern from the source inside the waveguide. 



Due to the resonator effect of the waveguide structure, 
the emission pattern does not show up on the scale of the 
contour plot in Fig. [Sj 
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FIG. 7: Contour plot of emission pattern, \G zz (r,r')\, of the 
system in Fig. [6] but for positions outside the photonic crys- 
tal. 



V. CONCLUSION 

We have described a procedure for solving the 
Lippmann-Schwinger equation for electromagnetic scat- 
tering in which the field along with the electric field 
Green's tensor is expanded in a basis of cylinder functions 
(so-called normal modes) inside each scatterer. Projec- 
tion of the electric field and the Green's tensor onto the 
normal modes is facilitated by the use of a number of ad- 
dition theorems to simplify the integral expressions and 
we have presented the method in general along with a 
thorough discussion of the evaluation of the scattering 
matrix elements, which may be helpful for practical im- 
plementation. 

The basis of normal modes ensures that all basis func- 
tions have the correct wave number. This, combined with 
the need for solving the system inside the scattering el- 
ements only, results in a relatively small linear equation 
system, as compared with other methods. Consequently, 
the method is fast and capable of handling large mate- 
rial systems such as photonic crystals. Furthermore, the 
use of a local cylinder function basis avoids the intro- 
duction of fictitious charges which may lead to instabil- 
ities for large refractive index contrasts in the case of 
TE polarization [l5[ and the integration scheme is free 
of staircasing errors along the boundaries. Due to the 
formulation in terms of the Green's tensor of the back- 
ground medium, there is no need for a calculation domain 
and the radiation condition is automatically satisfied, as 
are the boundary conditions (limited only by the numer- 
ical precision chosen). The accuracy of the method is 
thus limited only by the number of basis functions and 
the tolerance on the numerical integrals employed for the 
evaluation of the scattering matrix elements. We have in- 
troduced a measure of accuracy based on self-consistency 
and we have shown it to be of practical as well as principal 
importance. Once the matrix equation has been set up, 
it holds all information necessary to carry out scattering 
calculations on the geometry at the chosen frequency. It 
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can thus be stored and used for different choices of in- 
coming fields as well as for the calculation of the Green's 
tensor between different points r and r'. 

We have illustrated the method by two example prob- 
lems and we have shown an application of the method 
where we have calculated the zz component of the 
Green's tensor of a finite sized photonic crystal waveg- 
uide. Similar calculations will find application in the 
development of nanophotonic devices such as in the de- 
sign of junctions or cavities in photonic crystals or in- 
vestigation of the emission pattern from single photon 
sources. Using a similar procedure as the one described 
the method may be extended to three dimensional scat- 
tering geometries and although we have focused on ap- 
plications in micro- and nano photonic structures, we be- 
lieve the method may be of use in other areas of electro- 
magnetic scattering calculations as well. 
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FIG. 8: Sketch of relative coordinates as used in the expres- 
sion for Graf's addition theorem. 



We consider the cylindrical coordinates r = (r, tp) and 
r' = (r',<//), centered at two different positions O and 
O', respectively, where (L, 9) denotes the coordinates of 
O' with respect to O as shown in Fig. [8j Using this 
notation we express the Graf addition theorem as 



APPENDIX A: ADDITION THEOREMS FOR 
MULTIPOLE EXPANSIONS 

The expansion of the Lippmann-Schwinger equation 
used in this work, and especially the calculation of ma- 
trix elements, rely heavily on the use of cylinder func- 
tions. A number of addition theorems exist for cylinder 
functions which may simplify the calculations consider- 
ably. Of special interest in the present work is the Jacobi- 
Anger identity and the Graf addition theorem [3l| . Be- 
low we summarize the results in forms suitable for the 
present application. 

1. Jacobi- Anger identity 

For a plane wave, travelling at an angle 9 with respect 
to the x axis, the expansion in terms of cylinder functions 
is given by the Jacobi- Anger identity 

oc 

e t k or cos ( v -e) _ i n e- in8 J n (k r)e inv , (Al) 
n=— oc 

in which (r, ip) are cylindrical coordinates and J n {kor) is 
the Besscl function of order n. 



2. Graf's addition theorem 

The Graf addition theorem may be used to express 
the cylinder functions in one local coordinate system in 
terms of cylinder functions in a different local coordinate 
system. 



fl — — 00 

(A2) 

where Z n is a solution to Bessel's differential equation 
for integer n. If Z n = J„, the expansion is valid for all 
values of r', otherwise it is valid only for r' < L. 



APPENDIX B: CALCULATION OF MATRIX 
ELEMENTS 

In this appendix we evaluate the integral 

If= f G^(r,R)J,(k R R)(-ire-^dR, 

JP-8V 



which enters the expression for A a @ in Section III Bl 

For TM polarization, (a, ft) = [z, z), the angular inte- 
gration is nonzero only for \x = and the resulting inte- 
gral is well behaved, allowing for easy evaluation. For TE 
polarization, the integrand has a pole at the origin, so we 
rewrite this integral in a form more suitable for numeri- 
cal quadrature. Although the procedure is the same, the 
resulting integrals differ slightly depending on which of 
the elements of the Green's tensor we consider (see Ref. 
i29j for explicit expressions for the elements) . For I* x we 
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get 



TXX 



+ 





2i cos 



n 



sin 2 9H {k B R) + C °. S( " 2 ^ H±(k B (R) 



(20)1 



7T k%R 2 J 
■ r oc i-27, 

— lim 

4 5^0 



k B R 

J„(k R R)(-l)' J -e- il " e RdedR 



r'2it 2i cos(2#) 

/ i — L J u (knR)(— l)^e~ z ^ iid#di& a tes and Z x is a solution to Bessels differential equation 

./:: 7T k%R 2 



APPENDIX C: DERIVATIVES FOR GENERAL 
CYLINDER FUNCTIONS 



Below we summarize in polar coordinates, the 
double derivatives of the general cylinder function 
Z\(kr) exp{±iA</>}, in which (r, <fi) are cylindrical coordi- 



for integer A. 



ITT 
T 



(- 2<V -2 + <W - ^^ 2 )H {k B R) 



+ (<V -2 + <V,2A — ; — ^ h 



2i 



7Tfc|i? 2 



— lim , 

2 5-0 J 5 



(<5 M ,-2 + ^,2) 



Mk R R) 
k 2 B R 2 



RdR 



dx 2 



{Z x (kr)e iX *} = ^ [z x+2 (kr)e^ x+2 ^ 
+Z X - 2 {kr)e^- 2 ^ - 2Z x {kr)e iX *} , (CI) 



The first integral is now well behaved and may be di- 
rectly evaluated, whereas for the second integral we may 
use the identity 



lim 

5 



°Js 



k 2 B R 2 



RdR 



1 



In a similar way we rewrite the expressions for and 



If. v as follows: 



ITT 

TVV _ _ 

V 4 



-2 + ^ : 0^ M ,2)^o(fcs-R) 



1 

7T 
"4 

4fc 



2 (^-2 + <W0- 



-2 - S^ 2 )( 



lH 2 (k B R) 
2 k B R 



2i 



a t(\+2) V 



^ 2 {z x (kry^} = - k ^{z x+2 (kry 

+Z x _ 2 (kry( x - 2 ^ + 2Z x (kr)e lX ^ , (C2) 



dxdy 



irk 2 B R 2 



)Mk R R) RdR 



{Z x (kr)e iX *} = -i^ {z A+2 (fcr)e J ( A+2 )^ 

-Z A _ 2 (fcr)e l ( A - 2 ^}. (C3) 
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